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Abstract.  We  present  a  new  method  for  preprocessing  and  organizing  discrete  scalar  volume  data  of  any  dimen- 
sion on  external  storage.  We  describe  our  implementation  of  a  visual  navigation  system  using  our  method.  The 
techniques  have  important  applications  for  out-of-core  visualization  of  volume  data  sets  and  image  understanding. 
The  applications  include  extracting  isosurfaces  in  a  manner  that  helps  reduce  both  I/O  and  disk  seek  time,  a  priori 
topologically  correct  isosurface  simplification  (prior  to  extraction),  and  producing  a  visual  atlas  of  all  topologically 
distinct  objects  in  the  data  set.  The  preprocessing  algorithm  computes  regions  of  space  that  we  call  topological 
zone  components,  so  that  any  isosurface  component  (contour)  is  completely  contained  in  a  zone  component  and 
all  contours  contained  in  a  zone  component  are  topologically  equivalent.  The  algorithm  also  constructs  a  criticality 
tree  that  is  related  to  the  recently  studied  contour  tree.  However,  unlike  the  contour  tree,  the  zones  and  the  criti- 
cality tree  hierarchically  organize  the  data  set.  We  demonstrate  that  the  techniques  work  on  both  irregularly  and 
regularly  gridded  data,  and  can  be  extended  to  data  sets  with  nonunique  values,  by  the  mathematical  analysis  we 
call  Digital  Morse  Theory  (DMT),  so  that  perturbation  of  the  data  set  is  not  required.  We  present  the  results  of  our 
initial  experiments  with  three  dimensional  volume  data  (CT)  and  describe  future  extensions  of  our  DMT  organizing 
technology. 

Keywords:  volume  visualization,  isosurface  extraction,  contours,  Morse  theory,  data  organization,  data  compres- 
sion, digital  Morse  theory 


1.  Introduction 

Many  applications  produce  data  in  the  form  of  a  scalar 
function  defined  at  discrete  points  in  space  (called  vol- 
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ume  data).  Examples  include  X-ray  crystallography, 
Computed  Tomography  (CT),  Magnetic  Resonance 
Imaging  (MRI),  and  data  from  continuous  field  simu- 
lations such  as  computational  fluid  dynamics.  The  data 
is  visualized  using  either  volume  rendering  [1,  22,  23, 
26,  31]  or  isosurface  extraction  (see  [8,  9,  12,  21,  25]). 
In  volume  rendering  the  data  is  regarded  as  consisting 
of  semi-transparent  material  and  volume  primitives  are 
directly  projected  onto  the  screen  using  ray  casting. 
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Isosurface  based  methods  are  primarily  used  for  fil- 
tered volume  data.  I  n  i  sosurf  ace  based  methods  the  data 
is  segmented  by  thresholding,  that  is,  identifying  ob- 
jects defined  by  the  sets  of  points  in  space  where  a 
suitable  interpolating  function  /  for  the  filtered  data  is 
greater  than  or  equal  to  a  real  value  r,  called  the  iso- 
value  or  threshold.  The  topological  boundary  of  these 
objects  will  generally  correspond  to  the  level  sets  of 
/,  that  is,  the  points  p  where  f{p)  =  x.  The  data  is 
either  regularly  gridded,  that  is,  represents  a  function 
8  defined  at  points  of  Zn  on  the  vertices  of  a  hypercu- 
bic  decomposition  of  space,  or  it  is  irregularly  gridded 
and  S  is  typical  ly  defined  on  theverticesof  a  simplicial 
decomposition  of  Rn.  The  simplices  or  (hyper)cubes 
are  called  cells.  The  isosurface  extraction  problem  is 
to  construct  the  level  set  of  /,  called  the  isosurface,  for 
specific  isovaluer.  Individual  connected  components 
of  this  set  are  typically  called  the  contours,  even  when 
they  are  higher  dimensional  surfaces. 

M  any  researchers  have  developed  methods  of  orga- 
nizing the  data  so  that  the  "active  cells",  those  cells 
that  a  given  contour  intersects,  can  be  quickly  accessed 
from  disk  for  out-of-core  rendering.  The  three  main 
methods  include  geometric  organization  ([16,  34,  35], 
valuebased  organization  [4-7, 11, 24, 30],  and  recently 
the  topology  based  organization  using  augmented  con- 
tour trees  [2,  3,  32,  33]).  The  I/O  optimal  interval  trees 
of  [4-6]  seem  to  produce  the  least  disk  I/O  and  seek 
time,  and  input  the  minimal  number  of  blocks  in  active 
cell  acquisition.  Our  organization  is  topology  based. 

Topological  organization  of  thedata  uses  M  orseThe- 
ory.  M  orse  theory  studies  the  changes  in  the  topology 
of  the  level  sets  of  a  M  orse  function  /  as  the  parameter 
x  is  decreased.  These  changes  occur  at  the  values  of 
critical  points  of  /,  where  the  derivative  vanishes.  A 
M  orse  function  has  the  property  that  the  critical  points 
are  isolated  and  the  Hessian  is  nonsingular  at  each  of 
these  points.  Typically  the  data  readings  are  given  on 
a  tetrahedral  mesh  and  are  perturbed  to  insure  that  no 
two  val  ues  are  identical ,  and  a  si  mple  (e.g.  I  i near)  i nter- 
polant  /  is  selected,  to  insurethat  /  is  M  orse.  The  con- 
tour tree  tracks  the  topological  changes  of  individual 
contours  as  the  parameter  x  is  varied.  The  augmented 
contour  trees  of  [33]  are  used  to  produce  a  seed  cell 
set  and  active  cells  for  a  given  contour  are  acquired  by 
local  propagation  from  the  seed  cells. 

We  define  a  variant  of  the  contour  tree  and  use  it  to 
organizethedata  i  nto  topological  zones  and  topological 
zone  components.  We  demonstrate  that  the  zones  have 
important  properties  not  satisfied  by  the  augmented 


contour  tree,  making  them  useful  in  avariety  of  visual- 
ization problems.  For  example,  the  organization  aids  in 
active  cell  acquisition,  and  is  superior  to  the  seed  cell 
method. 

Prior  work  on  contour  trees  requires  that  regularly 
gridded  data  be  regridded  onto  tetrahedral  cells.  This 
has  two  disadvantages.  It  requires  increasing  the  data 
complexity  (e.g.  dividing  each  cubic  cell  into  5  tetra- 
hedra)  and  further  requires  an  explicit  representation 
of  the  tetrahedral  cells  and  their  adjacency,  which  in- 
creases storage  requirements  and  the  complexity  of 
dealing  with  the  data.  M  oreover,  the  data  readings  are 
then  perturbed  to  insurethat  no  two  readings  are  identi- 
cal.  I  n  this  way  one  i  nsures  that  thedata  can  be  extended 
to  a  M  orsefunction.  Our  mathematical  analysis,  which 
we  call  Digital  Morse  Theory,  demonstrates  that  this 
treatment  of  regularly  gridded  data  is  not  needed.  The 
contours  produced  by  the  most  popular  isosurface  ex- 
traction methods(seeforexample[28])  are  sufficiently 
well  behaved  so  that  the  topological  changes,  as  the 
isovalue  is  varied,  are  well  defined  and  easily  identi- 
fiable, without  recourse  to  M  orse  theory.  We  feel  that 
this  is  important  because  sampled  density  data  associ- 
ated with  a  physical  phenomenon  may  contain  regions 
of  constant  density,  and  regularly  gridded  data  does  not 
require  an  explicit  storage  of  cell  adjacency  informa- 
tion. We  give  a  combinatorial  characterization  of  the 
topology  changing  criticalitiesand  state  a  related  open 
problem. 

The  paper  is  organized  as  follows:  In  Section  2  we 
briefly  discuss  typical  methods  for  isosurface  extrac- 
tion, so  that  we  may  identify  the  properties  of  the 
surfaces  produced  by  these  methods.  We  also  give  a 
very  simple  algorithm  that  computes  contours  in  all 
dimensions  for  regularly  gridded  data.  In  Section  3 
we  give  our  mathematical  analysis.  We  formally  de- 
fine the  properties  satisfied  by  the  above  isosurface  ex- 
traction methods,  and  demonstrate  that  these  properti  es 
uniquely  determine  both  the  topology  and  the  topolog- 
ical changes  of  the  contours.  Alternatively,  one  can  use 
the  standard  tertrahedral  regridding  of  regularly  grid- 
ded data  and  perturbation  of  identical  readings  to  insure 
that  the  function  is  M  orse.  In  either  case,  the  zone  orga- 
nization isan  important  contribution,  and  the  rest  of  the 
paper  may  be  understood  independently  of  Section  3, 
with  the  usual  assumptions  of  a  M  orse  function. 

In  Section  4  give  the  formal  definition  of  the  topo- 
logical zones  and  prove  certain  of  their  properties.  In 
Section  5  we  show  how  to  compute  the  data  cells  that 
intersect  a  topological  zone  component  and  disussour 
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implementation  of  the  algorithm.  In  Section  6  we  de- 
scribe ourinitial  implementation  of  a  visual  navigation 
system  for  volume  data,  using  the  zone  organization. 
We  discuss  our  experimental  results  using  the  naviga- 
tion tool,  with  particular  attention  to  the  I/O  efficiency 
of  isosurface  extraction  using  our  technique.  We  also 
discuss  how  we  shall  scale  up  the  tool  to  handle  very 
large  data  sets.  In  Section  7  we  discuss  future  applica- 
tions of  our  technique,  including  topological ly  correct 
volume  data  simplification,  and  efficient  management 
of  level  of  detail.  Finally  in  the  appendix  we  give  the 
proofs  to  our  mathematical  assertions  of  Section  3. 


2.    Isosurface  Extraction 

The  isosurface  extraction  problem  is  often  stated  as 
follows.  A  scalar  volume  data  set  consists  of  tuples 
(x,g(x)},  where*  isa  point  in  3  (or  higher)  dimensional 
space,  and  g  is  a  scalar  function  defined  over  a  discrete 
set  of  these  points.  Given  an  isovalue.?,  the  isosurface 
extraction  problem  is  to  compute  the  (hyper)surface 
consisting  of  the  points  {x  :  f(x)  =  q},  where  /  is  a 
function  (interpolant)  which  extends  g  to  all  of  space. 
A  more  mathematically  precise  formulation  is  to  say 
that  we  are  computing  the  topological  boundary  of  {x  : 
fix)  >  q],  as  the  set  of  points  of  constant  valueg  may 
in  fact  include  a  3D  volume  in  the  degenerative  case, 
and  wemay  havemultiplesurfacecomponents,  which, 
by  an  abuse  of  terminology,  are  called  contours  in  the 
I  i  teratu  re.  A  f  u  rth  er  req  u  i  rement  i  s  th  at  th  e  co  n  to  u  rs  a  re 
oriented  manifolds. 

The  data  points  are  usually  organized  into  cells,  with 
the  data  readingsgiven  at  the  vertices  of  the  eel  Is.  Reg- 
ularly gridded  data  is  given  on  the  vertices  of  cubic 
cells  and  irregularly  gridded  data  is  typically  given  on 
the  vertices  of  a  tetrahedral  decomposition  of  space. 
The  algorithms  discussed  below  construct  the  isosur- 
face locally  within  each  cell. 


any  surface  intersects  the  cell,  and  then  moves  to  the 
adjacent  cell.  To  determine  surface  intersection  within 
a  cube,  the  algorithm  examines  all  its  vertices  and  as- 
signs a  High,  if  the  data  value  is  higher  than  or  equal 
to  the  isovalue,  (i.e.  these  vertices  are  considered  to  be 
inside  the  object)  or  Low,  if  the  data  value  falls  below 
the  i soval ue  ( i .e.  these  verti ces  are  outsi de  the  obj ect) . 

A  cubeedge  is  intersected  bythesurfaceif  andonly  if 
one  of  its  two  vertices  is  inside,  and  the  other  is  outside 
the  regi  on  of  space  bounded  the  surface  (the  i  sosurfaces 
are  assumed  to  be  oriented  manifolds).  The  location  of 
the  intersection  point  can  be  approximated  by  linear  in- 
terpolation. The  topology  of  the  surface  within  a  cube 
can  be  determined  by  the  pattern  of  its  vertex  values 
(H  igh  or  Low).  There  are  256  ways  to  color  8  vertices 
with  2  colors,  however,  by  taking  into  account  rota- 
tional and  complementary  symmetry  there  are  only  15 
topological^  distinct  patterns.  A  look-up  table  is  built 
that  includes  all  15  patterns  and  a  corresponding  sur- 
face approximation  for  each.  The  surface  approxima- 
tion uses  triangles  that  connect  the  intersection  points. 
All  the  distinct  patterns  of  surface  approximation  are 
enumerated  in  Fig.  1.  Each  High  is  illustrated  by  a  dark 
dot  at  a  vertex. 
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2.1.    Marching  Cubes 

A  classic  algorithm,  that  is  applied  when  the  volume 
data  is  given  on  cubic  cells,  was  proposed  by  Lorensen 
and  Cline  [25].  It  is  called  marching  cubes  and  it  con- 
structs a  polygonal  mesh  to  represent  a  3D  surface.  In 
order  to  generate  an  isosurface,  corresponding  to  a  user 
specified  isovalue,  the  algorithm  visits  all  cubic  cells 
created  from  eight  adjacent  pixels,  four  each  from  two 
adjacent  slices,  and  for  each  cell  determines  whether 


Figure  l.   Distinct  patterns  of  surface  intersection. 
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Figure  2.  An  example  illustrating  flaw  in  marching  cubes. 


Thegoal  is  obvi  ously  to  divide  each  cubi  c  cell  by  one 
or  several  surface  patches  (up  to  4  for  3  dimensional 
cells),  each  homeomorphic  to  a  unit  disk,  so  that  the 
induced  cell  regions  each  contain  only  connected  sets 
of  High  and  Low  vertices  (above or  below  theisovalue 
resp.)  Further  each  cell  face  contains  0,  1,  or  2  curve 
segments  that  divide  the  connected  H  igh  and  L  ow  ver- 
tices in  the  face.  The  5  possible  patterns  are  shown  in 
Fig.  9.  As  we  shall  see  this  goal  can  be  generalized  to 
higher  dimensional  cubic  cells. 


It  was  later  pointed  out  by  Nielson  and  Hamann  [27], 
that  the  original  marching  cubes  produced  tiling  errors 
between  cells,  and  he  proposed  an  asymptotic  decider 
to  remedy  this.  The  problem  occurs  in  what  Nielson 
terms  ambiguous  faces,  those  faces  with  4  separate  iso- 
surf ace  intersection  points  (which  weterm  4  hitfaces). 
Figure  2  shows  one  such  case  where  triangles  pro- 
duced by  marching  cubes  do  not  produce  a  continuous 
surface. 

For  instance,  if  a  cubic  cell  with  configuration  case 
6  shares  a  face  with  a  voxel  having  a  configuration 
3,  (which  is  the  complement  of  case  3),  triangles  pro- 
duced by  the  marching  cubes  algorithm  will  create  a 
hole  in  the  resultant  isosurface.  In  order  to  correct  this 
problem,  a  differenttri angulation  can  be  used.  Figure3 
depicts  two  possible  triangulations  that  will  result  in  a 
topological^  consistent  isosurface.  The  choice  of  tri- 
angulation  can  be  made  using  the  proposed  asymptotic 
decider  which  is  based  on  bilinear  interpolation  of  the 
value  of  an  interior  point  in  the  face. 


Figure  3.  Two  possible  triangulations  which  produce  correct  iso- 
surfaces. 


2.2.    Exploiting  Tetrahedral  Cells 

Irregularly  gridded  data  is  typically  organized  into 
tetrahedral  cells.  The  faces  of  the  tetrahedral  cells  can 
only  be  classified  (or  colored  )  in  three  distinct  ways, 
where  the  connectivity  among  the  vertices  can  be  re- 
solved without  any  ambiguity.  A  look-up  table  is  thus 
not  needed  to  find  the  approximated  surface.  The  three 
possible  cases,  described  in  Fig.  4,  are  the  following: 

-  only  vl  is  outside  the  object,  the  rest  are  inside. 

-  vl  and  v2  are  outside  the  object,  and  v3  and  v4  are 
inside. 

-  only  v4  is  inside  the  object,  the  rest  are  outside. 

Decomposing  cubic  cells  into  five  tetrahedral  cells 
[12]  can  lead  to  an  efficient  surface  extraction,  pri- 
marily because  of  the  fact  that  this  method  does  not 
suffer  from  the  well  known  inconsistency  discussed 
above.  However,  Zhou  et  al.  [37]  showed  that  tetra- 
hedral decomposition  and  linear  approximation  along 
the  introduced  diagonals  change  the  original  function 
and  may  lead  to  incorrect,  though  consistent  topology. 
They  proposed  trilinear  approximation  across  the  di- 
agonal of  the  cells  as  a  solution  to  this  problem,  rather 
than  applying  linear  approximation  along  all  edges. 
This  method  has  a  disadvantage  of  increasing  five- 
fold the  number  of  cells  that  have  to  be  processed 
individually. 

The  surface  construction  for  tetrahedral  cells  can  be 
generalized  to  higher  dimensions,  where  each  cell  isan 
^-dimensional  simplex.  In  each  case  the  vertices  above 
and  below  the  isovalue  in  a  cell  can  be  separated  by 
a  single  hypersurface  patch  that  is  homeomorphic  to 
an  n  -  1-dimensional  closed  ball.  For  cubic  cells  in 
n  dimensions  we  should  similarly  require  that  surface 
patches  within  the  cell  are  each  homeomorphic  to  the 
n  -  I  dimensional  unit  ball,  and  similarly  divide  the 
cell  into  regions  containing  connected  sets  of  High  or 
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Figure  4.   Surface  intersection  within  tetrahedral  cells. 


Low  vertices.  Thefol  lowing  simple  algorithm  achieves 
this  goal. 

2.3.    Spider  Web 

The  Spider  Web  algorithm  [9,  17,  21],  determines  a 
cell  intersecting  the  isosurfaceand  then  determines  hit 
points  along  the  cell  edges  using  linear  interpolation. 
In  the  current  implementation  of  the  algorithm  [10], 
the  pair  of  hit  points  on  a  voxel  face  with  2  hits  are 
considered  adjacent  and  a  decision  of  pair-wise  hit  ad- 
jacency on  a  voxel  face  containing  4  hits  is  made  (see 
Fig.  7)  using  a  strategy  similar  to  Nielson's  asymptotic 
decider. 

The  adjacency  defines  one  or  several  connected  sets 
of  hits  within  the  cell.  After  that,  the  centroid  point  of 
all  the  hits  within  a  connected  set  is  calculated.  This 
point,  called  the  articulation  point  (AP),  is  used  as  the 
cell  interior  triangle  vertex  for  all  of  the  subsequently 
constructed  triangles  within  that  set  of  hits.  Each  tri- 
angle consists  of  the  A P  and  a  pair  of  adjacent  hits 
on  a  cube  face  (see  Fig.  5).  Hits  on  a  cubic  cell  face 
are  shared  by  the  cell  that  shares  this  face,  and  thus 


Figure  5.   Constructing  isosurfaces  using  SpiderWeb. 


the  algorithm  can  continue  surface  construction  in  the 
adjoining  cells. 

Though  this  algorithm  has  the  advantage  of  not  re- 
quiring a  case  table,  it  produces  more  triangles  than 
marching  cubes.  On  the  other  hand,  it  is  completely 
paral lei izabl e,  as  each  cell  can  be  processed  indepen- 
dently. The  primary  cases  where  it  produces  more  tri- 
angles are  when  the  hit  set  consists  of  3  or  4  hits.  In  the 
former  case  one  can  remove  the  articulation  point  to 
produce  one  triangle.  When  there  are  4  hits,  marching 
cubes  produces  two  triangles.  H  owever  the  inclusion  of 
the  interior  A  P  and  the  subsequent  4  triangles  produces 
a  better  approximation  of  the  small  scale  curvature  of 
the  surface. 

Also,  unlike  the  original  marching  cubes,  it  is  eas- 
ily proven  correct.  We  sketch  the  proof  here  since  the 
ideas  are  used  in  the  proof  of  Theorem  1  (see  Fig.  6). 
To  prove  that  the  surfaces  produced  are  manifolds  one 
establishes  two  facts.  Each  traingle  edge  is  shared  by 
two  triangles,  and  the  set  of  triangles  incident  on  any 
triangle  vertex  are  homeomorphic  to  a  unit  disk.  The 
proof  uses  the  fact  that  each  hit  has  a  unique  adjacent 
hit  in  each  cell  face  on  which  it  occurs. 

The  triangle  edge  between  an  A  P  of  a  cell  C  and  a 
hit  h  is  shared  by  two  triangles:  one  triangle  for  each 
for  each  of  the  two  faces  of  C  that  share  this  hit.  Each 
triangle  consists  of  theAP,  h,  and  the  unique  adjacent 
hit  in  that  face.  The  triangle  edge  connecting  two  adja- 
cent hits  in  a  cell  face  is  shared  by  a  triangle  to  an  AP 
in  the  cell  that  shares  this  face. 

The  triangles  incident  on  an  AP  form  a  tri  angulation 
of  the  circle,  as  starting  from  any  hit  ft  in  the  connected 
set  and  following  hit  adjacency,  one  returns  to  h  and 
establishes  that  the  hits  form  a  cycle.  Similarly,  if  h  is 
any  hit  point,  h  is  incident  to  8  triangle  edges:  4  to  the 
unique  adjacent  hit  in  each  of  the  4  cell  faces  sharing 
the  cubic  cell  edge  on  which  h  occurs,  and  4  triangle 
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Figure  6.  (A)  The  triangle  edge  from  an  AP  to  a  hit  shared  by  2 
triangles.  (B)  The  edge  across  a  face  shared  by  2  triangles.  (C)  A  hit 
is  locally  a  tri angulation  of  the  circle. 


edges,  each  to  an  A  P  of  each  of  the 4  cubic  eel  Is  sharing 
the  eel  I  edge  on  which  h  occurs.  A  gai  n  these  8  tri  angles 
will  form  a  cycle. 

The  algorithm  can  be  extended  to  higher  dimen- 
sions. For  an  n  >  3  dimensional  cell,  recursively  build 
the  (hyper)  surface  patches  in  the  n  - 1  dimensional 
boundary  cells.  Each  of  these  patches  will  consist  of 
n  -  2  dimensional  simplexes  and  each  patch  will  be 
homeomorphic  to  an  n  -  2  dimensional  ball.  Now  for 
each  connected  set  of  hits  in  the  n  dimensional  cell 
select  an  A  P.  For  each  pair  of  adjacent  hits  in  the  set 
construct 2"-3  simplices,  each  consisting  of  theAP  and 
an  n  -  2  dimensional  simplex  that  contains  the  pair  in 
a  boundary  cell.  Each  such  simplex  will  consist  of  the 
pair  of  hits,  theAP,  and  n  -3  lower  dimensional  artic- 
ulation points. 

The  simplices  thus  formed  with  the  AP  for  the  set 
will  form  a  (hyper)  surface  patch  homeomorphic  to  an 
n  - 1  dimensional  unit  ball.  For  example,  in  a  4  dimen- 
sional cell,  the  lower  dimensional  patches  will  be  the 
normal  triangular  meshes  in  the  3  dimensional  bound- 
ary cubes.  For  a  connected  set  of  hits,  these  patches 
will  form  a  closed  surface  within  the  4  dimensional 
cell,  topological ly  equivalent  to  a  sphere  (by  the  proof 
of  correctness  of  the  original  Spider  Web).  Tetrahedra 
will  be  formed  from  an  interior  AP  to  each  triangle 
on  the  closed  surface,  forming  a  3  dimensional  surface 
patch  equivalent  to  a  ball.  This  patch  will  intersect  the 
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patches  in  neighboring  hypercubic  cells  by  completely 
sharing  a  2D  surfacepatch  in  a3D  boundary  cube,  since 
this  boundary  cube  will  be  completely  shared  by  two 
neighboring  4  dimensional  cells.  In  this  way  a  collec- 
tion of  3  dimensional  manifolds  is  constructed. 

Finally,  analyzing  this  behavior  of  Spider  Web  (and 
marching  cubes)  led  its  inventors  to  develop  the  new 
techniques  (Digital  M  orse  theory)  summarized  below 
[10, 18-20]. 

3.    Mathematical  Analysis:  Digital  Morse  Theory 
Explains  Why  Regridding  and  Perturbation 
are  Unnecessary 

The  basic  idea  of  the  analysis  is  as  follows.  We  charac- 
terize interpolants  that  produce  contours  topological^ 
equivalent  to  a  corrected  marching  cubes  (or  Spider 
Web  for  higher  dimensions)  or  the  standard  techniques 
for  tetrahedral  cells.  For  any  isovalue  one  labels  each 
data  reading  as  either  H  igh  (greater  than  or  equal  to  the 
isovalue)  or  Low  (less  than  the  isovalue).  The  topol- 
ogy of  the  isosurfaces  produced  by  the  aforementioned 
algorithms  is  uniquely  defined  by  this  binary  image, 
up  to  the  disambiguation  required  for  4  hit  faces.  The 
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behavior  is  such  that  the  surfaces  produced  induce  con- 
nected components  of  space.  We  call  the  objects  the 
components  that  are  i  nduced  by  (eel  I  edge  and  eel  I  face 
diagonal)  connected  sets  of  H  ighs  and  complementary 
objects  are  the  components  induced  by  connected  sets 
of  Lows  (here  connectivity  is  reminiscent  of  digital 
topology  [15]).  Thus  as  the  isovalue  is  varied  so  that 
it  passes  through  no  data  or  disambiguation  value  the 
topology  is  invariant.  As  the  isovalue  is  decreased  all 
that  can  happen  is  that  at  a  data  value,  one  or  several 
Lows  can  become  High,  or  at  a  disambiguation  value, 
two  diagonally  oppositeH  ighs  in  a  cell  facecan  become 
locally  connected.  Topological  changes  to  the  contours 
can  only  occur  in  the  following  ways. 

A  new  component  of  Highs  can  be  created  (anew  ob- 
jectisformed)  when oneorseveral  LowsbecomeHigh, 
or  a  component  of  Lows  can  vanish  (complementary 
object  destroyed).  Two  or  more  components  of  Highs 
can  be  merged  when  oneorseveral  LowsbecomeHigh 
or  connectivity  of  H  ighs  changes  at  a  disambiguation 
value,  resulting  in  separate  objects  merging.  Similarly, 
a  component  of  Lows  can  be  split,  resulting  in  com- 
plementary objects  splitting.  Finally  a  change  in  the 
connectivity  of  a  component  of  Highs  (or  Lows)  can 
induce  a  change  in  the  topological  type  of  one  or  more 
of  the  contours  bounding  the  object  (or  complementary 
object). 

The  important  point  is  that  one  can  infer  these 
changes  from  the  local  pattern  of  the  data  readings 
by  analyzing  the  behavior  of  the  isosurface  algo- 
rithms, without  any  recourse  to  M  orseTheory.  Wecan 
give  a  combinatorial  characterization  of  the  topology 
changes,  without  having  to  consider  the  analytic  prop- 
erties of  the  underlying  interpolating  function. 

3.1.  Definitions 

Definition  1.  We  will  assume  the  volume  data  set  is 
given  by  a  real-valued  function  s  defined  on  a  discrete 
set  of  points  V  in  R".  Regularly  gridded  data  is  given 
by  a  function  defined  on  Z",  where  we  form  unit  hy- 
percubes  in  the  natural  way.  Irregularly  gridded  data  is 
given  on  the  vertices  of  a  simplicial  decomposition  of 
Rn  and  we  assume  that  vertex  adjacency  information  is 
represented  in  some  form  in  external  storage.  We  term 
both  hypercubesand  simplices,  cells.  Without  any  loss 
of  generality,  we  shall  assume  S  is  non-negative  over 
the  entire  domain,  and  that  8  >  0  on  a  finite  sub-set  of 
points.  If  S  does  not  take  the  same  value  on  any  two 
points,  we  say  that  it  satisfies  data  uniqueness. 


For  simplicity,  we  will  first  describe  our  algorithm 
for  volume  data  that  satisfies  data  uniqueness.  Later 
we  will  extend  it  for  non-unique  data,  where  identi- 
cally valued  spatially  proximate  data  readings  imply 
isovolumes. 

Definition!.  A  continuous  real  valued  function  /  in- 
terpolates 8  if  it  extends  8  to  all  of  Rn  and  is  nonzero 
only  on  the  finite  subset  of  eel  Is  for  which  8  is  nonzero. 

We  denote  by  f~H>r),  the  set  of  p  e  Rn  such  that 
f(p)  >  t,  for  /  that  interpolates  8.  The  topological ly 
connected  components  of  this  set  are  called  objects. 
The  components  of  the  complement  of  this  set  are 
called  complementary  objects,  and  they  of  course  share 
common  boundary.  Similarly,  ^_1(>r)  denotes  the  set 
of  p  e  V  such  that  8{p)  >  r.  Clearly  if  /  interpo- 
lates 8  then  S-^^r)  c  f-H>r).  The  points  p  e  V 
whose  function  value  is  above  or  equal  to  (respectively 
below)  the  threshold  are  termed  Highs  (respectively 
Lows). 

Definition  3.  The  isosurface  construction  problem  is 
to  compute  the  contours  of  the  boundary  of  the  set 
/-1(>t),  for  specific  real  number  r.  We  use  the  nota- 
tion Bif"1^))  to  denote thetopological  boundary  of 
set  /_1(>r),  which  encloses  all  the  Highs  of  the  data 
region  for  isovalue  x. 

The  commonly  used  isosurface  extraction  meth- 
ods discussed  in  the  previous  section  construct 
#(/-1(>r))  with  the  simplest  topology,  consistent 
with  the  data.  These  methods  interpolate  a  single  con- 
tour intersection  point,  called  a  hit,  along  a  cell  edge, 
if  and  only  if  the  two  end  points  are  identified  as 
High  and  Low,  respectively.  Objects  are  defined  by 
connected  set  of  Highs.  For  irregularly  gridded  data, 
connectivity  between  the  Highs  is  defined  solely  by 
cell  edge  adjacency.  However  for  regularly  gridded 
data,  cell  edge  adjacency  is  not  sufficient,  as  observed 
above.  The  disambiguity  occurs  when  for  a  range  of 
thresholds  a  cube  face  F  contains  diagonally  oppo- 
site Highs  and  diagonally  opposite  Lows.  We  call  F 
a  4-hit  face,  because  there  are  hits  on  each  edge.  We 
will  assume  that  a  consistent  disambiguation  method  is 
chosen. 

Definition  4.  The  Disambiguation  value  c  with  re- 
spect to  a  function  /  is  the  maximum  threshold  for 
which  the  diagonally  opposite  Highs  in  F  are  locally 
connected  through  the  interior  of  a  cell  sharing  theface. 
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Note  that  once  the  data  points  are  path  connected 
they  will  remain  connected  for  all  threshold  less  than 
c.  We  associate  a  disambiguation  point  with  F ,  having 
thevaluec,  which  may  actually  be  in  thecube  interior, 
and  extend  S  to  this  point.  For  isovaluer  >  c  we  create 
a  pseudo-edge  connecting  the  two  diagonally  opposite 
Lows.  For  any  isovalue  r  <  c  a  pseudo-edge  connects 
the  two  diagonally  opposite  Highs.  Weextend  the  def- 
initions of  adjacency  and  connectivity  of  Highs  and 
Lows  to  include  the  pseudo-edges. 


3.2.    Properties  of  Admissible 
Interpolating  Functions 

Werestricttheclassof  interpolantsto  those  whosecon- 
tours  behave  in  a  manner  consistent  with  standard  iso- 
surface  extraction  methods  discussed  in  the  previous 
section.  The  following  are  the  properties  satisfied  by 
admissible  interpolants. 

Lett  bean  isovalue  not  in  the  range  of  the  extended 
8  (not  a  data  value  or  disambiguation  value). 

-  For  n  dimensional  data  the  contours  of  B  ( f  _1(>?)) 
form  a  finite  collection  of  disjoint,  compact  and  ori- 
ented n  -1  dimensional  manifolds  embedded  in  Rn. 
Thesecontoursdivide  Rn  intodisjointn  dimensional 
connected  components,  each  of  which  contains  ei- 
ther a  single  non-empty  connected  component  of 
Highs  or  Lows.  An  example  for  the  3  dimensional 
case  is  shown  in  Fig.  8,  where  the  space  is  divided 
between  disjoint  3  dimensional  connected  compo- 
nents of  Highs  and  Lows,  by  their  2  dimensional 
boundary  manifolds. 

-  The  intersection  of  the  contours  with  any  d- 
dimensional  cell  C  (d  <  n),  B(  f-1(>r))  n  C  ,  is  a 


Figure  8.   Connected  components  in  3D. 


finite  set  of  (d  -  1)  dimensional  components  called 
surface  patches.  These  surface  patches  are  bicon- 
tinuously  mappableto  a  (d  -  1)  dimensional  unit 
disk  (ball)  and  the  intersection  divides  C  into  dis- 
joint d  dimensional  contractible  components,  each 
of  which  contains  precisely  the verticesfrom  a  single 
non-empty  connected  component  of  Highs  or  Lows. 

Forr  in  therangeof  S  we  require  that  B  ( f  _1(>r))  = 
B  (Plx<r  f  _1(>x)).  This  insures  that  the  boundary  en- 
closes all  the  data  readings  precisely  equal  to  t. 

De/?nition5.  A  function  f  isadmissibleif  andonly  if 
i  t  i  nterpol  ates  Sand  sati  sfies  the  properti  es  I  i  sted  above. 

Property  1  merely  expresses  the  goal  of  isosurface 
construction  from  a  global  point  of  view,  and  Property  2 
expresses  the  type  of  local  non-degeneracy  assumed  by 
isosurfaceconstruction algorithms.  Forexample,  Fig.  9 
illustrates  all  possible  intersection  of  an  isosurfacewith 
a  cubic  cell  face  and  the  associated  division  of  the  cell 
face.  For  a  tetrahedral  cell  of  any  dimension  there  is 
always  single  patch  that  intersects  the  cell.  For  cubic 
cells,  thenumberof  distinctsurfacepatchescan  go  upto 
four.  Thissituation  occurs  when  thereexistsa  particular 
threshold  for  which  all  cell  faces  are  4-hit  faces  and  the 
given  isovalue  is  above  the  disambiguation  value  of 


Figure  9.   Possible  isosurface  intersections  with  cube  face. 
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all  faces.  Each  of  the  4  High  vertices  will  be  part  of  a 
disjoint  component  of  Highs. 

3.3.  Admissibility  Uniquely  Defines  Topology 

We  now  show  that  these  properties  are  sufficient  to 
uniquely  define  the  topology  of  irregularly  gridded 
data.  For  regularly  gridded  data,  the  disambiguation 
values  of  4-hitfacesarealso  required  to  uniquely  define 
thetopology  of  the  contours.  For  regularly  gridded  data 
the  theorem  i  s  real  ly  j  ust  a  tautol  ogy ;  i  sosurf aces  topo- 
logically  equivalent  to  SpiderWeb  or  M  arching  Cubes 
are  topological^  equivalent.  M  ore  precisely, 

Theorem  1.  If  f  and  f  are  admissible  and 
additionally,  for  regularly  gridded  data,  both  have  the 
same  disambiguation  values,  then  the  boundary  con- 
tours (manifolds) of  B ( f  _1(>t)) and  B(f ,_1(>r))are 
homeomorphic. 

Proof:  We  prove  this  theorem  by  a  series  of  lemmas 
(see  Appendix).  Here  is  an  outline  of  the  proof.  Wefirst 
observe  that  no  cell  facecan  havea  odd  number  of  hits. 
T  he  dual  of  data  reading  connectivity  is  the  hit  connec- 
tivity. We  identify  each  hit  by  the  cell  edge  on  which 
it  occurs.  Hits  are  adjacent  if  they  are  connected  by  a 
curve  segment  of  the  contour  in  a  cell  face.  Property  2 
ensures  that  a  pair  of  hits  on  a  cell  face  with  2  hits 
are  adjacent.  For  a  4-hit  face  the  hits  will  be  divided 
into  2  adjacent  pairs  (see  Fig.  7).  Since  both  f  and 
f '  are  admissible  and  have  the  same  disambiguation 
value  for  all  4-H it  faces,  they  will  have  the  same  hit 
set  and  the  same  connectivity  between  the  hits  for  any 
isovalue.  Let  M  and  M  '  be  manifolds  of  B(f-1(>T)) 
and  B  ( f  '~H>t}},  respectively.  We  form  a  graph  from 
the  surface  patches  of  a  manifold  where  the  nodes  of 
the  graph  are  patches  and  the  edges  connect  intersect- 
ing patches.  Since  intersecting  patches  must  share  a 
common  hit  and  the  patches  intersect  by  sharing  lower 
dimensional  patches,  itfollowsfrom  Property  2  that  we 
can  iteratively  construct  a  continuous  bijection  from  M 
to  M  '.  I  n  other  words,  one  can  show  that  the  graphs  G 
and  G ',  formed  from  M  and  M  ',  are  isomorphic.  □ 

3.4.  Critical  Values  for  Volume  Data  Satisfying 
Data  U niqueness 

For  the  rest  of  the  discussion  we  refer  to  components 
of  f_1(>r)  as  objects  and  components  of  f_1(<r)as 
complementary  objects.  Topology  changes  occur  when 


the  isovalue  becomes  equal  to  a  critical  value.  Let  us 
first  look  at  the  possible  critical  values  and  the  associ- 
ated changes  in  topology  before  arguing  how  the  ad- 
missibility uniquely  defines  topology  changes. 

L et us assumethat 8  satisfies data-uni queness.  A  crit- 
ical value  c  is  a  value  at  which  thetopology  of  thecon- 
toursof  f  change  as  the  threshold  isdecreased  through 
c.  Each  such  change  has  an  associated  critical  point  p 
for  which  f  (p)  =  c.  We  use  the  term  critical  point  in 
analogy  to  M  orse  theory,  but  we  are  not  actually  using 
any  differential  properties  of  the  function.  We  use  the 
admissibility  properties  of  the  function  to  identify  the 
topology  changes. 

Let  p  g  V  .  Let  N  be  the  set  of  points  in  V  adjacent 
to  p  by  cell  edgeorpseudo  edge.  Let  N  H  (respectively 
N  l)  bethe H ighs i n  N  (respectively  Lows)  with  respect 
to  the  value  of  p.  Compute  the  connected  components 
of  N  H  andNL  w  i  th  i  n  the  eel  I  s  that  share  p  butexcluding 
p  itself.  For  example,  two  Highs  are  connected  if  they 
are  connected  within  the  cells  that  share  p  by  a  path 
that  does  not  include  p.  Note  that  the  fourth  type  of 
criticality  only  occurs  in  cubic  cells. 

1.  I  f  N  =  NL,  meaning  p  has  a  higher  valuethan  all  its 
neighbors,  then  p  is  a  local  maximum  critical  point 
and  f  ( p)  is  a  local  maximum  critical  value.  A  new 
contour  emerges  at  each  of  these  critical  poi  nts  (see 
Fig.  10). 

2.  If  N  =  N  H ,  meaning  p  has  a  lower  value  than  all 
its  neighbors,  then  p  is  a  local  minimum  critical 
point  and  f  (p)  is  a  local  minimum  critical  value. 
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Figure  10.   Local  maximum  critical  point. 
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Figure  11.   Local  minimum  critical  point. 

An  existing  contour  vanishes  at  each  of  thesecritical 
points  (see  Fig.  11). 

3.  If  either  N  L  or  N  H  consists  of  more  than  one  con- 
nected componentthen  p  issaddlecritical  pointand 
f  (p)  is  a  saddlecritical  value.  Twoormorecontours 
may  merge  into  one,  contours  may  split  into  two  or 
more  pieces  or  a  change  to  the  topological  type  of 
one  of  the  contours  (genus  change  in  3D)  can  occur 
at  these  critical  points.  Figure  12  gives  an  illusta- 
tion  of  these  topological  changes  at  saddle  critical 
points. 

4.  Let  p  be  the  disambiguation  point  of  a  cell  face 
F  with  disambiguation  value  c.  If  the  Highs  of  F 
(respectively  Lows)  are  not  part  of  the  same  com- 
ponent within  the  cells  that  share  F  above  (respec- 
tively below)  the  disambiguation  value  then  p  is 
saddlecritical  pointand  c  is  a  saddle  value. 

Critical  points  of  type  1,  2  and  3  can  be  encountered 
both  forsimplicial  and  regularly  gridded  data,  whereas 
critical  points  of  type  4  can  occur  only  for  the  latter 
data  type.  Topology  changes  to  the  contours  can  only 
occur  at  critical  values. 


Theorem  2.  Let  f  be  admissible  and  S  satisfy  data 
uniqueness.  If  [ri,  r2]  is  a  critical  value  free  interval 
then  Btf-^Ti))  and  B ( f  ~1(>r2))  have  the  same 
topological  type  (homotopy  type).  The  converse  holds 
in  3-dimensions. 

See  the  A  ppendix  for  the  proof.  N  ote  that  we  have 
not  proved  the  converse  for  dimensions  n  >  3,  as  it  is 
beyond  the  scope  of  this  paper. 

3.5.   Critical  Isosets  When  Data  Does  Not 
Satisfy  U  niqueness 

If  the  data  does  not  satisfy  data  uniqueness  then  the 
properties  of  the  admissible  functions  imply  that  con- 
nected sets  of  data  points  (as  Highs)  of  identical  value 
c  will  be  enclosed  in  isovolumes  at  thresholds  c  -  e. 
As  the  threshold  is  decreased  through  c,  an  object 
can  merge  with  the  isovolume  containing  these  points. 
Complex  topology  changes  can  occur  at  the  valuec  of 
these  sets,  termed  i  sosets.  L  et  S  be  an  i  soset  w  i  th  val  ue  c 
and  N  ,  N  H  and  N L  bedefined asintheprevioussection. 
For  an  admissible  function  f ,  in  addition  to  the  point 
critical i ties  described  above  there  are  set  criticalities. 

1.  If  N  =  N  L ,  then  S  is  a  maximum  isoset  and  c  is  a 
maximum  critical  value.  A  new  contour  emerges  at 
each  of  thesecritical  isosets. 

2.  If  N  =  N  H ,  then  S  is  a  minimum  isoset  and  c  is  a 
minimum  critical  value.  An  existing  contour  van- 
ishes at  each  of  these  critical  isosets. 

3.  If  either  N  L  or  N  H  consists  of  more  than  one  con- 
nected component,  then  S  is  a  saddle  isoset,  andc  is 
a  saddle  critical  value.  Two  or  more  contours  may 
merge,  contours  may  spit  into  two  or  more  pieces, 
or  a  change  in  topological  type  may  occur  at  these 
critical  isosets. 


Figure  12.    Saddle  critical  point. 
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4.  If  N  H  and  NL  each  consist  of  a  single  nonempty 
componentthen  S  is  called  a  regular  isoset.  In  some 
cases  the  addition  of  the  isoset  to  an  object  (or  re- 
moval from  a  complementary  object)  may  cause  a 
change  in  the  topological  type  of  a  contour.  We  call 
this  a  critical  regular  isoset.  In  3  dimensions  this  can 
cause  a  genus  change  in  an  object  or  complimentary 
object  0  if,  when  it  merges  with  0  it  adds  or  de- 
stroys a  handle.  This  occurs  if  there  exists  a  loop 
through  N  L  or  N  H  and  a  loop  through  5  so  that  the 
loops  are  interlocked. 

Note  that  a  new  object  of  arbitrary  topological  com- 
plexity with  possibly  multiple  boundary  manifolds  is 
created  ata  maximum  isoset.  Oneorseveral  manifolds 
vanish  and  multiple  objects  can  merge  at  a  minimum 
isoset.  Extremely  complex  merges,  joins,  splits,  and 
Betti  number  changes  can  occur  at  saddle  isosets.  At 
a  regular  isoset,  in  3  dimensions  for  example,  a  genus 
change  can  occur.  The  conditions  are  that  the  structure 
of  the  isoset  itself,  when  surrounded  by  isosurfaces, 
is  topological^  complex  (not  simply  connected),  and 
additionally  it  is  joined  with  an  object  in  an  appropri- 
ate way.  For  example,  consider  a  donut-shaped  critical 
regular  isoset  with  flat  bottom  and  sides.  Suppose  this 
c  ri  ti  cal  set  merges  w  i  th  a  so  I  i  d  cu  be  o  bj  ect.  I  f  the  merge 
is  such  that  thedonut  rests  flat  on  top  of  the  cube,  then 
no  topological  change  occurs  because  boundary  of  the 
resultant  object  is  still  homeomorphic  to  a  sphere.  But 
if  the  merge  is  such  that  it  is  gl  ued  to  the  top  of  the  cube 
resting  on  its  side,  then  a  new  handle  will  be  formed 
causing  a  genus  change.  Both  of  these  cases  are  de- 
picted in  Fig.  13. 

We  leave  it  as  an  open  problem  to  develop  a  ef- 
ficient combinatorial  algorithm  for  recognizing  criti- 
cal isosets  in  all  dimensions.  Note  that  since  one  can 


Figure  13.   Non-critical  and  critical  regular  isoset. 


construct  a  si mplicial  complex  representing  the  isosur- 
faces one  can  in  fact  determine  algorithmically  whether 
a  regular  isoset  is  critical  (and  indeed  compute  both 
the  homotopy  groups  and  the  homology  groups  for  our 
contours!)  but  this  is  far  from  an  efficient  way  of  de- 
termining criticality.  This  is  not  a  practical  problem  as 
we  will  regard  regular  isosets  above  a  specific  size  as 
critical.  This  will  only  serve  to  possibly  increase  the 
number  of  zones  that  we  compute  (see  below)  but  will 
not  change  our  results.  The  important  point  is  that  topo- 
logical changes  can  only  occur  at  our  defined  critical 
values  and  ourinitial  experimentswithCT  dataindicate 
large  regular  isosets  are  rare. 

With  the  addition  of  the  critical  isosets  we  can  drop 
the  data  uniqueness  requirement  and  Theorem  2  can  be 
revised  as  follows  (seeAppendix  for  proof): 

Theorem  3.  Let  f  be  admissible,  if  [n,  r2]  is  a 
critical  value  free  interval  then  B(f_1(>Ti))  and 
B(  f  _1(>T2))  have thesame topological  type.Thecon- 
verse  holds  in  three  dimensions. 

H  ence  ami  ssibi  lity  uniquely  determines  the  topolog- 
ical changes  that  occur  as  the  isovalue  is  varied. 

4.    Topological  Zones  and  the  Criticality  Tree 
4.1.   Criticality  Tree 

The  criticality  tree  is  a  search  structure  that  hierarchi- 
cally organizes  the  data  cells  into  zones  based  on  value 
and  proximity,  and  traces  the  topological  evolution  of 
objects  in  the  data  as  the  isovalue  is  varied.  Each  node 
in  the  criticality  tree  is  a  criticality.  As  the  threshold  is 
decreased  from  the  value  of  a  critical  point  p,  if  q  isthe 
first  critical  point  encountered  in  the  same  object  that 
contains  p,  then  q  is  madethe  parent  of  p  (see  Fig.  14). 
M  ore  formally,  let  0  p(t)  be  the  object  containing  p  at 
threshold  r.  For  any  x  that  is  equal  to  a  data  value  define 
0  p(t)  to  be  the  closure  of  nc<T°  p'c)' 


Figure  14.  Topological  Zones  and  corresponding  Criticality  Tree. 
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De/inition 6.  Thenodesof thetreearethecriticalities. 
Let  p  be  a  critical  ity  and  letq  be  the  critical  ity  of  max- 
imum value  so  that  f  (q)  <  f  (p)  and  q  e  0  p(  f  (q)). 
Then  q  is  the  parent  of  p. 

If  the  data  does  not  satisfy  the  data  uniqueness,  q 
may  not  be  unique,  and  thus  the  object  containing  p 
may  meet  several  criticali  ties  at  the  same  time.  Wejust 
arbitrarily  order  the  parentage  of  these  to  break  ties. 

4.2.  Topological  Zones 

The  intuition  behind  the  topological  zones  is  quite  sim- 
ple. As  the  threshold  is  decreased  from  the  value  of 
a  critical  point  p,  the  object  containing  the  criticality 
growsand  deformswith  topological^  invariant  bound- 
ary until  it  encounters  another  criticality  q,  which  in 
some  manner  changes  the  topology  of  the  contours 
bounding  theobject.  Thetopological  zoneof  criticality 
p  is  the  volume  swept  by  the  topological  boundary  of 
this  object  until  it  encounters  q.  Each  connected  com- 
ponent of  this  zone  is  the  volume  swept  by  an  individual 
boundary  manifold  (contour)  of  the  object.  Figure  14 
gives  a  visual  representation  of  the  zones. 

We  define  the  zones  by  a  set  difference  so  that  no 
assumption  that  the  contours  vary  continuously  is  re- 
quired, because  if  the  data  contains  an  isoset  the  con- 
tours will  not  vary  continuously.  Thetopological  zone 
for  criticality  p  is  defined  as  follows: 

De/inition  7.  Let  p  be  a  critical  point  with  critical 
value  xp  and  parent  q.  Thetopological  zone  of  p,  de- 
noted ^(p),  is  the  set  difference  Op(f(q))  —  Op(f(p)). 
The  topological  zone  components  are  the  connected 
components  of  £(p). 

4.3.  Properties  of  the  Topological  Zones 

Theorem  4.  Let  p  be  a  criticality  with  parent  q  with 
a  =  f  (p)  and  b  =  f  (q). 

1.  For  each  x  satisfying  a  >  r  >  b,  B(0  p(r))  is  con- 
tained in  f(p)  and  each  component  M  of  B(Op(t)) 
is  contained  in  a  single  component  of  f(p).  This 
means  that  for  any  isovalue  between  the  value  of 
a  critical  point  and  that  of  its  parent,  the  contour 
is  completely  contained  within  the  zone  for  that 
criticality. 

2.  For  an  pair[ri,T2],  if  a  >  x\  >  r2  >  b,  B(Op(ti)) 
and  B  (0  p(r2))  are  topol ogi cal I y  equivalent. 


3.  For  all  x,  if  0  is  any  component  of  f  _1(>r)  not 
containing  p,  then  B(O)  and  f(p)  are  disjoint. 

4.  For  all  x,  such  that  r  >  a  and  r  <  b,  B(f_1(>r)) 
and  £(p)  are  disjoint. 

Proof:  A  point  r  e  £(  p)  if  and  only  if  f(q)<  f(r)  < 
f  (p)  and  r  e  0  p(  f  (q)).  For  any  r  such  that  f  (p)  > 
r  >  f  (q),  each  point  r  e  B(0  p(r))  satisfies  f  (r)  =  x 
by  continuity  of  f.  Also,  because  each  such  r  e 
0 p( f (q )),  B(Op(t))  is  contained  in  £(p)  by  mono- 
tonicity  of  f  _1(>r).  The  connectivity  of  the  manifolds 
insures  that  it  is  contained  in  a  single  component  of 
£(p),  establishing  Claim  1.  By  definition,  £(p)  con- 
tains no  critical  point  with  a  value  between  f  ( p)  and 
f  (q),  which  establishes  the  second  claim.  Claim  3  fol- 
lowsfromtheconnectivity  of  theO  p(r)  and  thefact  that 
boundary  manifolds  are  pair-wise  disjoint  by  Property 
1  of  admissiblef  unctions.  Claim  4  follows  immediately 
from  the  proof  of  Claim  1.  □ 

The  criticality  tree  and  zones  organize  the  data  set 
in  a  hierarchical  fashion,  in  the  following  sense: 

Corollary  1.  The  union  of  all  zones  in  the  subtree 
rooted  at  node  p  is  completely  contained  in  the  class 
of  objects  corresponding  to  p.  That  is,  for  any  isovalue 
r  in  the  range  of  f  ( p),  0  p(r)  consists  of  a  portion  of 
£(p)  and  the  union  of  all  zones  in  the  subtree  rooted 
at  p. 

4.4.   Criticality  Tree  vs.  Contour  Tree 

There  are  some  essential  difference  between  the  criti- 
cality tree  and  the  contour  tree  [2,  3,  32,  33].  The  two 
search  structures  are  contrasted  and  compared  below. 

-  The  contour  tree  traces  the  evolution  of  individual 
contours  whilethecriticality  tree  traces  the  evo I  uti  on 
of  the  objects  the  contours  bound. 

-The  contour  tree  does  not  record  a  change  in  the  ho- 
motopy  type  (genus  in  3D)  of  an  individual  contour, 
while  the  criticality  tree  records  both  a  change  in  the 
topological  type  of  any  contour  bounding  an  indi- 
vidual object,  as  well  as  a  change  in  the  number  of 
contours  bounding  an  individual  object. 

-  Leaf  nodes  in  the  contour  tree  can  be  maxima  or 
minima  where  contours  are  created  or  destroyed, 
as  the  threshold  is  decreased,  while  leaf  nodes  in 
the  criticality  tree  are  all  maxima,  where  objects  are 
created. 
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F  igure  15.   Level  set  of  f  as  the  parameter  x  decreases. 

-  In  a  contour  tree  a  node  with  multiple  children  can 
occur  when  two  or  more  distinct  objects  merge  or 
when  a  contour  of  an  object  splits.  In  a  critical  ity  tree 
multiple  children  only  occur  when  objects  merge. 
Objects  cannot  split  as  the  isovalue  is  decreased, 
though  complementary  objects  can  split,  resulting 
in  an  increase  in  the  number  of  contours  bounding 
an  object. 

-As  one  descends  the  criticality  tree  from  the  root 
the  isovalue  is  strictly  increasing.  The  root  of  the 
criticality  tree  is  the  single  object  that  encloses  the 
entire  data  set.  Following  a  path  from  the  root  to 
a  leaf  traces  the  evolution  of  a  single  object  as  the 
isovalue  is  increased.  All  the  objects  in  the  subtree 
of  a  node  p  are  completely  contained  in  the  root 
object  corresponding  to  p.  There  is  no  particular 
order  when  following  a  path  in  the  contour  tree. 

-The  seed  cell  method  [33],  which  uses  the  contour 
tree  for  finding  the  active  eel  Is,  extracts  each  contour 
by  searching  for  a  seed  cell  known  to  have  inter- 
sected the  isosurface  in  the  corresponding  super-arc 
first  and  then  locally  propagates  from  the  seed  cell. 
However,  for  an  extremely  large  data  set  this  local 
propagation  can  incur  excessive  disk  seek  time.  The 
zone  organization  method  uses  a  more  sophisticated 
method  to  reduce  I/O  operations.  Since  a  zone  com- 
ponent file  contains  all  the  active  eel  Is  for  a  contour, 
they  can  be  extracted  by  reading  thezonecomponent 
sequentially,  reducing  disk  seek  time  to  finding  the 
start  of  the  component. 

Let  us  illustrate  the  difference  between  the  con- 
tour tree  and  the  criticality  tree  through  an  example. 
Figure  15  shows  the  level  set  of  a  function  f  as  the 
parameter  t  is  decreased,  and  Fig.  16  shows  the  corre- 
sponding contour  tree  and  the  criticality  tree.  Initially 
there  are  four  objects  created  at  maxima  7  through 
10.  These  merge  at  saddles  5  and  6.  Then  the  con- 
tours bounding  each  of  the  objects  undergo  two  genus 
changes,  changing  from  sphere  to  torus  to  sphere  again 
at  saddles  5'  and  5"  ( at  saddles  6'  and  6"  respectively). 


Figure  16.   Contour  and  Criticality  Tree. 


These  changes  are  not  reflected  in  the  contour  tree. 
Then  the  two  objects  merge  at  saddle  4,  the  boundary 
of  the  this  object  becomes  toroidal  at  saddle  4'  and  a 
inner  boundary  is  created  at  saddle  3  by  a  split.  At  this 
point  the  object  comes  to  end  osea  hoi  low.  Note  that  the 
contour  tree  has  two  edges  coming  outof  node3,  which 
represents  the  outer  and  inner  boundary  manifolds  of 
the  object,  respectively.  Finally  the  bubble  is  closed  at 
minima  2  and  this  node  represents  the  single  object 
containing  the  entire  data  set. 

5.    Zone  Organization  Algorithm 

5.1.   Computing  Zones 

The  input  to  the  preprocessing  algorithm  is  a  volume 
data  set.  For  regularly  gridded  data  the  data  set  is  orga- 
nized as  a  stack  of  two-dimensional  slices.  Each  slice 
contains  a  two-dimensional  array  of  real  numbers.  The 
output  is  a  criticality  tree  and  a  direct  access  file  con- 
taining a  sorted  list  of  the  cells  (voxels)  intersecting 
each  zone.  The  criticality  tree  indexes  the  zone  file, 
with  each  node  containing  pointers  to  the  components 
for  the  corresponding  zone.  The  pre-processing  of  the 
volume  data  to  compute  the  tree  and  the  cells  that 
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intersect  each  zone  component  is  performed  in  three 
phases,  which  are  described  below.  The  zone  compu- 
tation algorithm  proceeds  by  assigning  a  label  to  each 
data  poi  nt,  w  hich  corresponds  to  the  zone  that  i  ncl  udes 
that  particular  point. 

-  Phase  1:  In  the  first  phase  the  data  is  scanned  and 
critical  points  (and  critical  isosets)  are  identified  by 
computing  the  connected  components  of  H  ighs  and 
Lows  in  the  neighborhood  of  each  point  and  identify 
isosets  by  depth-first  search.  One  assigns  the  label 
l(p)  =  p,  where  p  is  a  critical  point  (or  representa- 
tive point  of  a  critical  isoset).  A  zone  is  created  for 
each  critical  ity  p. The  critical  points  are  then  placed 
on  a  max-heap  ordered  by  value,  with  priority  given 
to  critical ities  to  break  ties. 

-  Phase  2:  One  then  iterates  the  following  step  until 
the  heap  is  empty.  Remove  the  maximally  valued 
point  r  from  the  heap  with  label  l(r)  =  p.  Point 
r  can  either  be  a  critical  or  non-critical  point.  We 
elaborate  the  actions  taken  in  both  the  cases  below. 

•  Non-Critical  Point:  Fora non-critical  pointr  ±  p, 
examine  the  neighbors  of  r.  If  aneighborq  of  r  is 
unlabel  led  then  assign  a  tentative  label  to  q  from 
the  label  of  r  (i.e.  I(q)  =  l(r)  =  p)  and  place  q 
on  the  heap.  The  label  of  r  is  now  finalized  and  it 
is  placed  in  the  zone  for  critical  ity  p. 

•  Critical  Point:  Foracritical  pointr  =  p,  examine 
theneighborsof  r .  If  there  is  any  unlabel  led  neigh- 
bor assign  the  tentative  label  as  before.  If  there  is  a 
neighbor  q  that  has  a  finalized  label  different  from 
the  label  of  r,  or  in  other  words  I (q)  =  a,  such 
that  a  ^  p,  then  make  the  node  that  corresponds 
to  the  criticality  a  the  parent  of  the  node  that  cor- 
responds to  the  criticality  r.  The  flooding  of  the 
zone  for  criticality  r  is  now  over.  The  points  that 
are  still  waiting  in  the  heap  with  a  tentative  label 
p  are  relabeled  tentatively  by  a,  and  each  one  of 
them  are  also  placed  in  thefilefor  criticality  r  (see 
discussion  below  for  efficient  relabeling).  Critical 
disambiguation  points  are  regarded  as  adjacent  to 
all  the  four  vertices  of  the  ambiguous  cell  face. 

-  Phase  3:  One  creates  a  single  direct  access  file  for 
all  the  zones.  One  next  iterates  through  the  original 
zone  files  dividing  each  zone  file  into  its  connected 
components  using  a  standard  depth  first  search.  Af- 
ter the  zone  components  are  computed  for  a  zone, 
a  cell  list  is  computed  for  each  component  and  out- 
put to  the  direct  access  file.  Each  zone  component 
stores  the  cell  information  (e.g.  coordinate,  values 


at  all  vertices,  local  maximum  and  local  minimum) 
of  al  I  the  eel  I  s  that  i  ntersect  the  component.  T  he  eel  I 
information  is  sorted  by  thelocal  maximum  first  fol- 
lowed by  the  local  minimum  values  in  order  to  able 
to  scan  fewer  cells  during  isosurface  extraction. 

Observe  that  the  algorithm  labels  the  entire  data  set 
by  highest  value  first  order.  That  ensures  that  a  zone  is 
expanded  as  much  as  possible,  by  expanding  the  zone 
boundary  in  the  shallowest  descent  direction,  before 
it  is  terminated.  Also,  at  any  point  in  time,  all  points 
greater  than  or  equal  to  a  particular  r  have  been  la- 
beled. T  hese  are  precisely  that  points  that  comprise  the 
components  of  f_1(>t).  Whenever  a  zone  labeling 
for  a  point  p  terminates,  the  points  remaining  in  the 
heap  with  a  tentative  label  p  are  the  ones  that  lie  in  the 
cells  that  contain  the  boundary  of  this  zone,  and  the 
data  points  with  the  finalized  label  of  p  are  precisely 
the  points  contained  in  f(p).  Cells  that  straddle  one  or 
several  zone  boundaries  are  placed  in  each  such  zone. 

5.2.   Implementing  the  Zone  Organization  Algorithm 

Our  primary  focus  is  limiting  the  I/O  and  maintaining 
as  little  in-core  data  as  possible  throughout  the  execu- 
tion. The  most  expensive  step  is  the  relabeling  of  the 
boundaries  of  a  zone  from  phase  2.  Instead  of  immedi- 
ately relabeling  the  points  we  use  a  lazy  approach  that 
results  in  correct  but  efficient  relabeling.  The  points 
that  have  to  be  relabeled  are  left  on  the  heap.  When 
a  point  is  removed  from  the  heap,  we  trace  the  par- 
ent pointer  of  the  node  that  corresponds  to  its  current 
label,  up  to  the  root  node  in  the  partially  built  critical- 
ity tree,  and  give  a  finalized  label  that  corresponds  to 
the  root  node.  While  following  the  parent  pointers,  the 
point  is  also  placed  in  all  the  zone  files  corresponding 
to  the  nodes  on  the  path.  The  cell  in  which  this  point 
lies  straddles  the  boundary  of  each  such  zone.  The  al- 
gorithm thus  runs  in  time  0  (kn  log  n),  where  n  is  the 
number  of  data  cells  and  k  is  the  longest  path  traversed 
in  thetreeduring  relabeling  (our  experiments  with  CT 
data  have  shown  this  number  to  be  small  in  practice). 

During  phase  1  the  maximum  number  of  slides  that 
are  read  into  main  memory  at  any  point  in  time  is 
restricted  to  three.  Suppose  that  currently  (i  -l),i, 
(i  + 1)  are  read  into  memory.  All  critical  points  in  the 
ith  slide  are  detected  and  then  the  next  slide,  (i  +2), 
is  read  into  the  same  space  as  (i  - 1).  If  the  data  in- 
cludes an  isoset  then  additional  data  needs  to  be  read 
from  arbitrary  slides.  In  that  case  only  the  necessary 
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data  points  are  brought  into  memory,  not  the  whole 
slide.  Execution  of  any  phase  can  be  carried  out  in- 
dependently given  that  the  previous  phases  have  been 
completed.  We  have  implemented  the  zone  preprocess- 
ing on  a  14  processor  Sun  Enterprise  U I traS pare,  and 
it  has  been  parallelized.  The  implementation  is  done  in 
an  out-of-core  fashion,  so  that  it  may  be  scaled  up  to 
larger  volume  data  sets.  During  Phase  1,  only  portions 
of  at  most  3  si  ices  are  kept  in  memory  at  one  time.  Dur- 
ing Phase  2  only  the  data  readings  on  the  heap  are  kept 
in  memory.  The  amount  of  in-core  data  is  proportionate 
to  the  contour  size.  We  also  have  a  mode  where  we  la- 
bel only  one  zone  at  a  time  for  further  space  efficiency. 
If  we  find  that  when  we  go  to  gigabyte  scale  datasets, 
the  memory  requirements  for  the  heap  grow  too  large 
(which  we  don't  anticipate,  see  below)  we  can  use  an 
external  memory  implementation  of  the  heap.  Phase  3 
requires  marking  of  visited  vertices  during  depth  first 
search  and  marking  of  cells  during  cell  list  construc- 
tion. We  presently  use  an  in-core  boolean  array.  How- 
ever further  space  efficiency  (at  the  cost  of  time)  could 
be  achieved  by  using  an  external  array  or  sparse  matrix 
technology  for  this  marking.  The  sorting  of  the  cells 
within  each  zone  is  done  by  an  efficient  sort.  The  sort- 
ing time  dominates  the  phase  and  thus  the  total  time 
for  Phase  3  is  bounded  by  0  (kn  log  kn),  where  k  is 
(as  above)  the  maximum  number  of  zones  any  point 
appears  in. 

Since  we  do  not  have  a  exact  procedure  to  iden- 
tify critical  regular  sets  we  assume  all  regular  isosets 
above  an  arbitrary  size  to  be  critical.  However  in  that 
case,  some  zones  may  be  divided  by  regular  isosets  that 
are  not  actually  critical ities.  This  will  not  change  the 
properties  of  the  zones  that  we  have  defined,  but  may 
segment  the  zone  unnecessarily  and  make  it  smaller. 
Again  we  have  found  large  regular  isosets  to  be  rare, 
even  in  unfiltered  CT  data. 


6.    Applications  of  the  Zone  Organization: 
A  Visual  Navigation  Tool 

We  have  implemented  a  visual  navigation  tool  using  the 
criticality  tree  as  a  search  data  structure,  with  pointers 
to  the  zones  and  the  zone  components  on  disk.  It  al- 
lows us  to  visualize  volume  data,  varying  both  the  iso- 
val  ue  and  the  spati  al  vi ew poi  nt  i  nteracti  vel y  i  n  near  real 
time.  This  tool  is  used  to  construct  the  contours  at  any 
given  isovalue  x.  The  tool  first  uses  the  criticality  tree 
to  find  the  active  zone  components.  The  program  scans 
through  the  cell  list  for  each  component.  Because  the 


cell  lists  are  sorted,  as  soon  as  a  cell  with  a  max  value 
max(C )  <  t  is  encountered,  the  search  for  active  cells 
is  terminated  for  that  zone  component.  The  SpiderWeb 
algorithm  contructs  the  surface  patches  within  each  ac- 
tive cell  independently,  and  thus  the  construction  can 
be  parallelized. 

The  zone  organization  allows  us  to  produce  a  visual 
atlas  of  all  topological ly  distinct  objects  in  the  data  set 
together  with  the  range  of  isovalues  that  reveals  each 
obj  ect.  A  common  cri  tici  sm  of  the  i  sosurf  ace  extracti  on 
methods  is  that  they  obscure  the  underlying  structure 
of  the  scalar  field.  The  zones  reveal  it.  The  visualiza- 
tion tool  also  allows  us  to  increment  or  decrement  the 
isovalue  by  any  value  entered  by  the  user.  The  tool 
reconstructs  the  isosurface  from  the  loaded  zone,  and 
only  loads  a  new  zone  when  critical  values  are  passed. 
The  tool  can  also  be  used  to  animate  the  evolution  of 
th e  o  bj  ects  as  the  i  so va I  u e  i  s  dec reased  f  ro  m  th e  g  I  o ba I 
maximum  to  global  minimum  by  a  small  increment. 
Figure  17  shows  few  snapshots  of  such  an  animation 
done  on  a  si  mulated  data  set  that  has  a  si  mi  lar  structure 


Figure  17.  Snapshotof  thelevel  setof  afunction  f  as  the  parameter 
r  is  decreased. 
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Table  1.   Experimental  results. 


Exp.  1 

Exp.  2 

Exp.  3 

Data  description 

Simulated  data 

Filtered  medical  data 

Filtered  medical  data 

Size  (Row  x  Col  x  Slide) 

40  x  40  x  40 

20  x  20  x  20 

100  x  100  x  100 

Zone  size  (as  %) 

13.95 

16.27 

3.09 

Percentage  of  active  cells 

28.83 

11.66 

2.21 

Percentage  of  cells  loaded 

42.6 

15.49 

2.94 

Load  utility  (%) 

67.67 

75.28 

75.17 

Travel  utility  (%) 

82.07 

94.87 

96.08 

to  thedata  set  shown  in  Fig.  15.Thesystem  also  allows 
the  user  to  rotate  the  objects  generated  at  any  isovalue, 
thus  the  user  can  interactively  search  for  the  isovalue 
that  reveal  s  a  structure  of  i  nterest.  T  he  tool  ki  1 1  oads  and 
visits  only  a  small  percentage  of  the  total  cells.  This  is 
done  on  a  low-end  SGI  workstation. 


6.1.    Experimental  Results 

We  have  experimented  on  simulated  data,  and  both  fil- 
tered and  unfiltered  CT  data  taken  from  the  Visible  Hu- 
man data  set  and  the  results  are  summarised  i  n  Table  1. 
Unfiltered  data  produces  many  more  zones.  We  need 
to  make  modifications  to  reduce  the  number  of  zones 
for  unfiltered  data.  If  thezone  is  small  enough  (e.g.  the 
size  of  thezone  is  less  than  or  equal  to  one  disk  block) 
we  can  read  this  single  zone  into  memory  using  one 
disk  seek. 

Experiment  1  was  conducted  on  simulated  data  that 
modeled  the  objects  in  the  screen  snapshots  shown  in 
Fig.  17.  In  this  way  we  could  verify  by  hand  the  cor- 
rectness of  the  zone  computation.  Experiment  2  was 
conducted  on  a  small  portion  (20  x  20  x  20)  of  un- 
filtered medical  data,  followed  by  a  filtered  (standard 
Gaussian)  version  of  the  same  data  set.  We  iterated 
through  the  isovalue  range  by  a  variety  of  small  incre- 
ments, and  verified  that  al  I  active  eel  Is  for  each  contour 
(for  each  isovalue)  were  contained  in  a  zone  compo- 
nent. Experiments  wasconducted  with  a  larger  portion 
of  thefiltered  CT  data  set  (100  x  100  x  100).  Wefound 
that  the  unfiltered  data  produced  too  many  zones.  We 
did  not  count  the  I/O  used  to  read  the  criticality  tree 
file  into  core,  as  this  is  done  once  at  the  start  of  the 
program.  Additionally,  as  we  iterated  through  the  iso- 
values,  in  many  cases  no  I/O  was  performed,  as  the 
necessary  zone  was  already  loaded.  Nevertheless,  for 
experimental  purposes  we  counted  the  I/O  that  would 


be  performed  in  loading  thezone.  Forthesepreliminary 
experiments  we  recorded  the  average  zone  component 
size  as  a  percentage  of  the  data  set  size,  the  average 
number  active  cells  (over  the  range  of  isovalues)  as  a 
percentage  of  the  data  set,  and  the  average  percentage 
of  cells  that  were  actually  brought  into  memory.  Note 
that  for  the  simulated  data  the  average  number  of  ac- 
tive cells  exceeds  the  zone  size,  as  at  most  thresholds 
several  zones  were  needed  to  rendertheobjects(azone 
contains  a  single  object). 

Finally  the  two  most  important  statistics,  in  ana- 
lyzing the  performance  of  our  active  cell  acquisition, 
are  the  load  utility  and  the  travel  utility.  The  load 
utility  records  the  average  percentage  of  active  cells 
among  those  brought  into  memory.  In  other  words,  it 
records  the  percentage  of  loaded  cells  that  were  actu- 
ally utilized  to  construct  the  isosurf  ace.  The  travel  util- 
ity records  the  percentage  of  active  cells  among  those 
actually  examined.  This  is  perhaps  the  most  important 
measure  of  efficiency.  Since  the  zones  are  sorted,  we 
stopped  examining  cells  in  a  zone  once  the  next  cell 
was  out  of  range  for  the  isovalue.  Since  our  goal  was 
to  reduce  disk  seek,  we  read  the  entire  zone  compo- 
nent. However  the  travel  utility  indicates  the  percent- 
age of  active  eel  Is  among  those  actual  ly  exami  ned;  we 
can  load  into  memory  only  these  examined  cells  if  we 
choose.  As  we  can  see  both  utility  factors  turned  out 
to  be  quite  high.  It  will  be  interesting  to  compare  these 
statistics  with  the  other  search  structures,  particularly 
the  interval  trees  of  [4]. 

The  average  percentage  of  in-core  data  is  around 
20%  for  small  portions  of  data,  which  is  quite  reason- 
able, because  on  average  15%  of  the  total  cells  are 
active  cells  that  we  actually  need  to  produce  the  iso- 
surf ace.  When  we  increase  the  sample  size  to  100  by 
100  by  100  (for  filtered  data)  the  percentage  drops  far 
below  10  For  filtered  CT  data  the  percentage  of  the 
cells  actually  examined  that  were  active  (part  of  the 
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Figure  18.  Rearrangement  of  the  Criticality  Tree  when  a  zone  is 
tagged  deprecated. 


surface)  averaged  between  94.87  and  96.08  percent. 
As  we  increase  the  data  set  size  we  expect  that  the 
zone  size  will  grow  at  about  n2/3.  When  we  increased 
the  data  set  from  about  8000  cells  to  about  1,000,000 
we  found  that  the  percentage  of  data  readings  in  a  zone 
dropped  more  than  five-fold. 

For  both  filtered  and  unfiltered  medical  data,  many  of 
zones  included  fewer  than  0.1  percent  of  the  total  data 
points,  Such  zones  do  not  really  reflect  any  significant 
change  inthetopology  andean  be  culled  without  signif- 
icant loss  of  information.  A  good  heuristic  for  culling 
would  be  the  number  of  new  data  points  included  in  a 
zone,  as  opposed  to  the  total  zone  size,  which  can  in- 
clude many  points  on  the  boundary  of  several  zones.  It 
w  i  1 1  be  i  nteresti  ng  to  experi  ment  w  i  th  usi  ng  the  zone  al  - 
gorithm  to  reduce  noise,  as  opposed  to  Gaussian  filter- 
ing, which  may  obscure  essential  features  (see  Fig.  18). 

6.2.   Scaling  U  p  the  Visualization  Tool 

The  topology  of  the  boundary  manifold  of  a  zone 
component  remains  invariant  for  any  isovalue  that  lies 
withi  n  the  range  of  that  component.  This  means  that  as 
I  ong  as  we  retai  n  the  cri  ti  cal  i  ti  es  we  can  coal  esce  eel  I  s 
within  a  zone  component  up  to  the  limits  given  by  the 
zone  structure  without  changing  the  topology  of  the  re- 
sulting contours.  In  this  manner  wecan  reducethe  tiling 
complexity  either  prior  to  or  during  extraction,  based  on 
the  users'  demand.  Thus,  wecan  have  finely  rendered 
surfaces  in  the  area  of  interest,  and  coarsely  rendered 
surfaces  in  other  areas.  This  is  different  from  simpli- 
fication after  construction,  see  for  example  [14,  29]. 
The  prospect  of  dynamic  multi resolution  (as  opposed 
to  preprocessed  multiresolution  [13]),  for  cubic  cells 
is  particularly  exciting  for  the  following  reason.  Once 
one  has  read  all  the  cells  from  a  zone  into  memory, 
refining  the  resolution  can  be  done  without  further  I/O 


Figure  19.  Coarsely  rendered  or  finely  rendered  without  changing 
the  topology. 


penalty.  For  example  one  can  form  larger  cubic  cells 
(say  64  by  64  by  64)  from  the  original  cells.  To  refine 
the  rendering  detail  recursively,  one  only  needs  to  in- 
clude an  interior  data  reading  in  the  larger  eel  I,  dividing 
it  into  8  smaller  cells.  Since  this  reading  is  already  in 
memory  no  I/O  is  required.  For  static  multiresolution, 
one  can  store  each  zone  component  at  several  scales. 
If  the  zone  component  is  small  enough,  all  the  scales 
can  be  brought  into  memory  initially.  If  not,  each  scale 
can  be  read  as  needed.  Figure  19  shows  an  example  of 
the  different  scales.  We  have  not  yet  employed  these 
multiresolution  methods,  but  they  will  be  essential  in 
producing  a  very  large  data  set  visualization  tool  with 
real  time  capabilities. 

7.   Future  Applications  and  Conclusion 

Itwill  be  very  interesting  to  explore  the  potential  useof 
topological  zones  in  applications  such  as  multi-modal 
image  registration  (see  [13,  14,  36]  for  good  surveys 
of  these  problems).  Note  that  the  criticality  tree  is  in- 
variant with  respect  to  rotation,  scaling,  translation  and 
uniform  value  translation  of  the  data  set.  Two  images 
with  same  contents  rotated  at  different  angle,  scaled 
differently  or  translated  will  result  in  similar  or  par- 
tially si  mi  lar  criticality  trees.  Thecriticality  trees  can  be 
compared  to  verify  whether  or  not  the  images  have  the 
same  content.  The  image  registration  problem  can  then 
be  reduced  to  registration  of  relevant  portions  of  the 
two  trees,  and  the  corresponding  critical  points. 

The  zone  segmentation  provides  a  very  natural  and 
efficient  way  of  organizing  the  data  set  si  nee  the  organi- 
zation is  based  on  thecontents  of  the  data  set.  The  whole 
procedure  of  the  preprocessing  is  automatic  and  does 
not  need  any  user  intervention.  The  zone  segmentation 
preprocessing  results  in  efficient  I/O  operation  which 
provides  a  great  starting  point  for  overcoming  the  bot- 
tleneck that  remains  one  of  the  unresolved  challenges 
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of  visualization  and  animation  of  very  large  data  sets. 
Sincethezonesegmentation  is  hierarchical  and  thecrit- 
icality  tree  traces  the  evolution  of  objects  we  believe  it 
can  be  used  as  a  modelling  tool  for  volume  data.  The 
nodes  of  the  subtree  rooted  at  any  node  p  of  the  crit- 
icality  tree  contain  all  the  cells  in  the  interior  of  the 
object  associated  with  p.  Thus  thetree  itself  providesa 
topological^  complete  model  of  all  threshold  defined 
objects  in  the  data  set. 

We  have  demonstrated  the  feasibility,  correctness 
and  utility  of  the  zone  organization  in  visualizing  vol- 
ume data  sets.  We  believe  that  we  have  provided  a 
basis  for  exploiting  this  powerful  technique  that  will 
ultimately  result  in  the  ability  to  visually  navigate  ex- 
tremely large  volume  data  sets  in  real  time  without 
special  purpose  hardware. 

Appendix 

Proof  of  Theorem  1: 

Lemma  1.  Let  e  be  a  cell  edge.  Then  e  has  a  single 
intersection  point  (hit)  with  B  ( f  _1(>t))  if  and  only  if 
the  incident  end  points  are  H  igh  and  Low,  respectively. 

Proof:   This  follows  immediately  from  Property  2. 

□ 

Lemma  2.  Each  2  dimensional  hypercube  face  F 
(square)  has  0,  2,  or  4  hit  points.  Each  2  dimensional 
simplex  face  F  (triangle)  has  0  or  2  hits. 

Proof:  This  follows  from  a  simple  parity  argument 
or  by  examining  the  4  possible  arrangements  of  Highs 
and  Lows  on  the  4  vertices  of  a  cube  face.  1,  2  or  3 
edge  adjacent  H  ighs  give  2  hits  and  a  pair  of  diagonally 
opposite  Highs  gives  4  hits  (see  Fig.  9).  □ 

D  e/inition  8.  Two  hits  are  adjacent  in  a  face  F  if  they 
are  connected  by  a  curve  segment  of  B  ( f _1(>r))  n  F . 
Connectivity  of  sets  of  hits  is  defined  by  the  transi- 
tive closure  of  adjacency.  Connected  components  are 
defined  accordingly. 

Lemma  3.  Each  hit  point  h  of  face  F  is  connected 
by  a  curve  segment  of  B  ( f  _1(>r))  n  F  to  a  unique  hit 
point  h'  in  F ,  that  is  adjacent  to  h. 

Proof:  The  result  follows  from  the  previous  lemma, 
the  disambiguation  rule  and  Property  2  (see  Fig.  9). 

□ 


De/inition  9.  We  uniquely  identify  a  hit  point  by  the 
cell  edge  on  which  it  occurs. 

Lemma  4.  B(f  _1(>t))  and  B(f'-1(>r))  have  the 
same  hit  sets  and  hit  connectivity  if  f  and  f '  are  both 
admissible. 

Proof:  Changing  the  interpolation  function  does  not 
change  whether  a  vertex  is  High  or  Low.  For  regularly 
griddeddata,  f  and  f '  both  havethesamedisambigua- 
tion  values  resulting  in  the  same  hit  connectivity.  □ 

De/inition  10.  Two  surface  patches  are  adjacent  if 
they  intersect. 

Lemma  5.  Two  d  -  1  dimensional  surface  patches 
of  M  ,  a  e  C  and  a'  e  C  are  adjacent  if  and  only  if 
their  intersection  is  a  set  of  d  -  2  dimensional  surface 
patches  of  M  in  the  d  -  1  dimensional  boundary  cell 
shared  by  C  and  C '.  F  urther,  if  a  and  a'  are  adjacent 
then  they  share  at  least  one  hit  point  in  each  component 
of  their  intersection. 

Proof:  T  he  fact  that  they  must  intersect  i  n  some  lower 
dimensional  patch  follows  from  the  fact  that  the  d  -1 
dimensional  boundary  cell  is  completely  shared  by  C 
and  C  and  Property  2  of  the  admissible  function.  The 
fact  that  they  share  a  hit  point  in  each  component  of 
thei  r  i  ntersecti  on  i  s  due  to  the  fact  that  thei  r  i  ntersecti  on 
must  occur  on  at  I  east  one  eel  I  edge,  and  thus  they  share 
the  hit  along  this  cell  edge.  □ 

Lemma  6.  Each  connected  component  of  hits  in 
cell  C  is  a  member  of  a  unique  surface  patch  ae 
B(  f  _1(>r))  nC  .Conversely  the  hits  contained  in  each 
surface  patch  a  e  B( f  _1(>r))  n  C ,  form  a  unique 
connected  component  of  hits  in  C  . 

Proof:  By  definition  of  hit  adjacency  and  connec- 
tivity, each  pair  of  hits  in  a  connected  component 
of  tine  hits  within  C  is  connected  by  a  path  within 
B  ( f  _1(>t))  n  C  .  By  the  path  connectivity  of  the  sur- 
face patches,  thi  s  path  must  be  a  part  of  the  same  patch. 
Thus  the  connected  set  of  hits  all  lie  in  a  single  patch 
a  of  C . 

For  the  other  direction,  we  must  show  that  the  hit 
poi  nts  of  a  form  a  connected  set  withi  n  C  .  The  proof  is 
by  induction  on  the  dimension  d.  The  result  is  clearly 
trueford  =  1  and  d  =  2.  For  inductive  hypothesis,  we 
assume  that  for  all  j  <  d  -1,  the  hits  of  each  patch  of  a 
j  dimensional  cell  C  is  a  connected  set  withi n  C .  Leth 
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and  h'  beany  two  hit  points  of  o.  Si  nee  a  ishomeomor- 
phic  to  the  unit  ball,  its  boundary  p  is  connected  and 
containedinthed-ldimensional  boundary  eel  Is  of  C  . 
Then  there  is  a  path  n  e  p  from  h i  to  h 2  such  that  this 
path  7T  passes  through  a  sequence,  01,  o2, . . . ,  om,  of 
adjacent  d  -2  dimensional  patches  in  thed  -1  dimen- 
sional boundary  cells  of  C  .  Let  hj  beany  hit  point  in  oj 
and  hj+i  be  a  hit  point  in  oj+i,  where  1  <  i  <  m.  We 
claim  that  hj  and  hi+i  are  part  of  the  same  component 
of  hits  of  C. 

To  provetheclaim  we  first  observe  that  si  nee  a,  and 
ffi+i  are  adjacent,  by  Lemma  5  they  share  at  least  one 
hit  point.  By  the  induction  hypothesis,  the  hits  of  Oj 
and  CTj+i  are  each  connected  individually  within  their 
respective  boundary  cells.  Thus  the  union  of  the  hits 
of  01  and  CTj+i  must  be  part  of  the  same  connected  set 
of  hits  in  C .  This  also  implies  that  hi  and  hj+i  are 
part  of  the  same  connected  component  of  hits  in  C  , 
thus  establishing  the  claim.  Any  two  hit  points  in  o  are 
part  of  the  same  connected  set.  In  other  words,  the  hit 
points  contained  in  any  surface  patch  o  form  a  unique 
connected  component  of  hits  within  the  cell  C  .  □ 

Lemma  7.  Each  connected  component  of  hits  is 
contained  in  a  unique  component  manifold  M  of 
B  ( f  _1(>r))  and  conversely  the  hits  contained  in  each 
such  manifold  M  form  a  unique  connected  component 
of  hits. 

Proof:  A  connected componentof  hitsiscontained  in 
a  single  manifold  M  by  the  definition  of  hit  adjacency 
and  hit  connectivity.  This  manifold  is  unique,  since  by 
Property  1  the  manifolds  are  pair-wise  disjoint. 

For  the  other  direction,  suppose  that  M  contains  two 
hit  points  hi  and  h2.  There  is  clearly  a  path  jt  in  M  from 
hi  and  h2.  This  path  passes  through  some  sequence  of 
patches  of  M  ,  01,05, ... ,  om  of  cells  C 1,  C  2,  •  ■  ■ ,  Cm, 
where  hi  e  01  and  h2  e  am.  For  each  1  <  i  <  m,  a\ 
contains  a  unique  connected  component  of  hits  within 
C  i  by  L  emma  6.  F  or  each  such  i ,  Oj  and  o\  +1  i  ntersect 
in  a  shared  boundary  cell  of  Ci  and  Ci+i  and  share  at 
least  one  hit  point  by  Lemma  5.  Thus  the  hit  points  of 
oj  and  oj+i  are  part  of  the  same  connected  component 
of  hits,  which  means  hi  and  h2  are  part  of  the  same 
connected  set  of  hits  as  well.  n 

From  the  previous  lemma  we  get 

Lemma  8.  Iff  and  f' are  both  admissible,  B(f_1 
(>t))  and  B(  f'-1(>T))  consistof  the  same  number  of 
components. 


Let  H  be  a  connected  component  of  hits.  Let  M 
and  M  '  be  the  components  of  B(f  _1(>t))  and  B(  f /_1 
(>r))  containing  H  . 

Lemma  9.  Let  M  and  M  '  consistof  the  patches  £  = 
01,  o2, . . . ,  ok  and  £'  =  a{,  o2, . . . ,  a'y,  respectively. 
There  exists  a  one-to-one  mapping  f  of  £  onto  £'  so 
that  for  all  1  <  i  <  j  <  k,  oj  is  adjacent  to  oj ,  if 
and  only  if  f(a\)  isalso  adjacentto  ^(oj).  In  addition 
their  respective  intersections  have  the  same  number  of 
components. 

Proof:  The  mapping  f  associates  two  patches  if  the 
share  the  same  hit  set.  The  mapping  is  one-to-one  and 
onto  by  Lemma  6.  If  a\  is  adjacent  to  a\ ,  then  from  the 
above  lemmas,  they  share  at  least  one  hit  in  common. 
Suppose  they  share  a  hit  on  cell  edge  e.  Then  f(o\ )  and 
Vr(oj)  also  share  a  common  hit  one  cell  edge  e  and  are 
adjacent.  □ 

In  other  words  if  we  represent  the  patches  and  their 
adjacencies  by  graphs  in  the  natural  way,  then  the 
graphs  representing  M  and  M  '  are  isomorphic.  Using 
the  above  fact  we  can  prove 

Lemma  10.    M  and  M  '  are  homeomorphic. 

Proof:  We  can  construct  a  global  homeomorphism 
from  M  and  M  '  piece-wise  via  their  patches.  From 
Property  2,  both  patches  oj  and  f(<j\)  are  bicontinu- 
ously  mappableto  a  unit  disk.  From  Lemma  9  we  know 
thatoj  is  adjacent  to  oj  if  and  only  if  f(o\)  isalso  ad- 
jacent to  f[o\).  Also  by  Lemma  5  and  Property  2  we 
know  that  the  i  ntersecti  ons  of  the  correspondi  ng  pai  r  of 
the  adjacent  patches  arebicontinuously  mappable.  We 
can  thus  bi continuously  mapoj  to  ^(01 )  and  extend  the 
mapping  for  each  adjacent  patch  recursively.  Thus  we 
can  define  the  homeomorphism  iteratively  on  each  oj 
and  extend  it  to  define  a  global  homeomorphism  from 
M  to  M  ',  proving  the  lemma.  □ 

Theorem  1  now  follows  immediately  from 
Lemmas  7  and  10. 

Proof  of  Theorem  2:  Now  the  contours  that  we  con- 
sider are  all  manifolds  and  each  can  begiven  by  asim- 
plicial  complex.  For  our  simple  case  weak  homotopic 
equivalence  implies  homotopic  equivalence:  it  is  suf- 
ficient to  know  that  the  homotopy  groups  of  two  con- 
tours are  isomorphic  to  give  topological  equivalence. 
What  this  means  is  that  we  need  consider  two  types 
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of  topological  changes  to  our  contours:  change  in  the 
number  of  contours  or  change  in  one  of  the  homotopy 
groups  (and  corresponding  Betti  number)  of  a  contour. 
Observe  that  objects  and  sets  of  Highs  are  monotonic 
in  the  sense  that,  as  the  threshold  is  decreased,  High 
data  readings  remain  High,  and  points  that  were  part 
of  the  object,  remain  a  part  of  the  object  throughout. 
The  Lows  becomes  High  after  certain  time  and  new 
points  are  added  to  the  objects.  F  urther,  from  the  proof 
of  Theorem  1  we  know  that  no  topological  changes  (in 
the  strong  sense  of  homeomorphism)  can  occur  as  the 
isovalue  is  varied  in  an  interval  that  does  not  contain  a 
data  or  disambiguation  value. 

Topological  changes  to  B  ( f  _1(>t))  can  thus  occur 
in  several  ways  and  of  course,  several  changes  can  oc- 
cur at  the  same  time.  First  of  all  a  new  object  can  be 
created  as  the  threshold  r  is  decreased  through  a  value 
c,  and  hence  a  new  contour  is  created.  It  is  straightfor- 
ward to  see  from  Property  1  of  admissible  functions 
that  this  can  only  happen  if  a  new  component  of  Highs 
iscreated,  and  thus  c  isalocal  maximum  critical  value. 
I  n  this  case  there  is  no  obj  ect  i  n  some  regi  on  of  spacef  or 
all  thresholds  t  >  c,  buta  new  object  exists  in  thesame 
regionforr  <  c.  By  admissibility  and  data  uniqueness, 
this  new  object  must  contain  a  single  data  point  p  and 
must  be  compri  sed  of  a  si  ngl  e  connected  component  of 
Highs.  Thus  p  is  a  local  maximum  critical  point.  Sim- 
ilarly a  contour  can  vanish  at  c,  and  by  monotonicity 
this  can  happen  only  if  c  is  a  local  minimum  critical 
value,  where  the  associated  point  p  comprises  a  single 
connected  component  of  Lows  for  x  =  c  +  e. 

Two  or  more  separate  contours  can  merge  at  valuec. 
F  rom  Property  land  monotonicity  this  meansthat  there 
exist  two  separate  components  of  Highs  at  threshold 
t  =  c  +  €,  that  are  merged  into  one  component  of 
Highs atr  =  c— e. Thiscan  only  happen  if  c  isa saddle 
critical  value  of  type  3  and  a  point  p  becomesHigh  ate 
and  unites  the  two  components,  or  c  is  a  critical  saddle 
point  of  type  4  and  the  change  of  the  connectivity  of 
H  ighs  in  a  4-hitface  unites  the  two  separate  component 
of  Highs. 

One  or  several  contours  can  split  at  c.  From  admis- 
sibility Property  1  and  monotonicity  this  will  happen 
only  when  a  complementary  object  (the  region  contain- 
ing a  component  of  Lows)  is  separated  by  a  change  in 
connectivity  at  c.  This  can  happen  only  of  c  is  a  saddle 
critical  value  of  type  3  or  4. 

In  3  dimensions  a  contour  can  change  in  genus  at  c, 
when  a  class  of  incontractible  loops  is  created  or  de- 
stroyed. In  higher  dimensions  either  a  homotopy  class 


of  incontractible  k-spheres  (k  <  n)  is  created  or  de- 
stroyed. Monotonicity,  data  uniqueness  and  admissi- 
bility imply  that  an  increase  in  a  Betti  number  of  a 
contour  can  only  occur  in  the  following  way.  As  the 
threshold  is  decreased  through  a  valuec,  separate  por- 
tions of  the  sameobject  that  existin  a  neighborhood  of 
a  point  p  at  threshold  c  +  e,  are  merged  at  threshold 
c  -  e,  so  that  a  new  class  of  incontractible  loops  or 
k-spheres  is  created.  The  fact  that  the  object  is  locally 
separated  in  the  cells  containing  p  at  c  +  e  implies 
from  admissibility  that  there  must  beat  least  two  com- 
ponents of  Highsinthese  cells,  thus  insuring  that  p  isa 
saddle  critical  point.  Symmetrically,  from  admissibil- 
ity a  decrease  in  a  Betti  number  can  occur  only  when  a 
complementary  object  is  locally  separated.  Again  this 
means  that  a  component  of  Lows  is  locally  separated 
at  a  value  c.  From  admissibility  and  data  uniqueness 
this  means  that  either  a  single  Low  becomes  High  or 
there  is  a  change  in  connectivity  in  a  4-hitface.  In  the 
former  case  the  data  poi  nt  that  changed  wi  1 1  be  a  saddl  e 
critical  point,  and  in  the  latter  case  the  disambiguation 
point  of  the  4-hitface  will  be  a  saddle. 

We  now  demonstrate  that  topological  changes  do 
occur  at  critical  values  in  3  dimensions.  By  admissi- 
bility, a  new  object  is  created  at  a  local  maximum  and 
a  complementary  object  vanishes  at  local  minimum. 
M  ultiple  changes  can  occur  at  saddles,  including  any 
combination  of  contour  merges,  splits,  and  handles  cre- 
ated or  destroyed.  We  argue  that  at  least  one  of  these 
changes  occur  at  a  saddle  value.  Suppose  that  saddle 
point  p,  with  value  c  is  part  of  the  same  global  com- 
ponent of  Lows  at  threshold  c  +  e  and  at  threshold 
c  -  €,  when  p  becomes  H  igh,  the  component  is  sepa- 
rated. By  Property  1  a  boundary  contour  has  split  ate. 
Similarly  if  two  globally  separate  components  of  Highs 
are  merged  at  c,  at  least  one  pair  of  the  corresponding 
boundary  contours  has  merged. 

To  demonstrate  that  changes  to  topological  type 
(genus  of  a  contour)  occur  at  saddle  values  when  ob- 
jects are  not  merged  (respectively  complementary  ob- 
jects are  not  split),  note  that  the  eel  Is  N  sharing  saddle 
point  p  (with  value  c)  are  homeomorphic  to  unit  ball 
B,  with  p  at  the  center  and  the  other  data  points  on 
the  spherical  boundary  of  the  ball.  Suppose  that  there 
are  at  least  two  components  of  H  ighs  of  N  H ,  N  i  and 
N2,  on  the  boundary  of  B  at  thresholds  above  c  (the 
argument  for  NL  is  symmetric).  By  assumption  both 
Ni  and  N2  are  part  of  the  same  obj  ect  0  .  We  claim  that 
one  can  find  a  loop  L  on  the  boundary  of  B  (on  the  eel  I 
faces)  that  separates  N 1  and  N2  on  the  boundary  of  B 
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so  that  L  is  exterior  to  0  .  The  definitions  of  connec- 
tivity of  Highs  and  Lows  insure  that  we  can  construct 
L  as  a  loop  consisting  of  cell  edges  and  cell  face  diag- 
onals connecting  adjacent  Lows.  0  will  intersect  the 
boundary  of  B  in  at  least  two  components  Oi  and  O2 
containing  Ni  and  N2  respectively.  By  connectivity  of 

0  and  its  boundary  contours,  one  can  find  a  path  it 
in  the  boundary  of  0  from  Oi  to  O2  so  that  n  never 

1  ntersects  the  i  nteri  or  of  the  eel  I  s  B ,  and  so  that  at  c  -  e, 
it  can  be  extended  to  a  loop  L '  by  a  path  from  0 1  to  0  2 
through  the  interior  of  B.  Clearly  L'  is  incontractible 
as  it  will  be  interlocked  with  the  above  loop  L.  □ 

Proof  of  Theorem  3:  As  in  the  proof  of  Theorem  2, 
from  admissibility  a  new  object  can  only  be  created  at 
the  value  of  a  maximal  isoset  or  point  and  a  comple- 
mentary object  can  only  be  destroyed  at  the  value  of  a 
minimal  point  or  isoset.  A  change  in  a  Betti  number  of 
a  contour  can  now  occur  in  two  ways.  In  the  first  way 
two  locally  separate  portions  of  an  object  can  merge  at 
a  value  (resp.  a  complementary  object  is  locally  sep- 
arated at  a  value)  in  such  a  way  that  a  new  homotopy 
classof  incontractible  loops  or  k-spheres,  k  <  n  is  cre- 
ated (resp.  destroyed)  on  its  boundary.  From  admissi- 
bility this  means  that  a  globally  connected  component 
of  Highs  is  merged  at  a  point  or  isoset  (respectively 
component  of  Lows  is  locally  separated)  and  this  im- 
plies that  the  isoset  or  point  is  a  saddle  criticality  as 
we  have  defined.  The  second  way  that  the  change  can 
occur  is  caused  when  an  isoset  is  merged  with  another 
object  (or  removed  from  a  complementary  object)  so 
that  the  structure  of  the  isoset  itself  adds  a  homotopy 
class  of  incontractible  loops  (or  spheres)  to  the  bound- 
ary of  the  object  (or  destroys  a  class  of  incontractible 
loops  or  spheres  in  the  boundary  of  the  complemen- 
tary object).  The  isoset  in  this  case  is  a  critical  regular 
isoset. 

We  now  argue  that  topological  changes  do  indeed 
occur  at  critical  values  in  3  dimensions.  A  change  in 
topological  type  occurs,  by  definition,  at  the  value  of 
a  critical  regular  isoset.  A  new  component  of  Highs 
is  created  at  the  value  of  a  maximal  isoset  and  hence, 
from  admissibility,  a  new  object  is  created  with  at  least 
one  new  boundary  contour.  Similarly,  a  component  of 
Lowsvanishesatthevalueof  a  minimal  isoset  and  thus 
a  complementary  object  and  at  least  one  boundary  con- 
tour is  destroyed.  If  two  separate  components  of  H  ighs 
are  merged  at  a  saddle  isoset  then  the  corresponding 
objects  are  merged  with  a  corresponding  change  in  the 
number  of  contours.  If  a  component  of  Lows  is  split 


by  a  saddle  isoset  then  a  new  complementary  object  is 
created  and  hence  the  number  of  contours  changes. 

We  now  consider  a  saddle  in  which  neither  of 
these  two  things  happens.  So  assume  without  loss  of 
generality  that  N  H  consists  of  at  least  two  connected 
sets  of  Highs  (the  case  of  Lows  is  symmetric)  that  are 
globally  part  of  the  same  object  0  above  the  critical 
threshold  c.  As  in  the  previous  proof,  let  N  be  the  eel  Is 
containing  the  saddle  set.  Divide  N  into  its  connected 
components  by  cell  edge  adjacency.  If  the  components 
of  Highs,  Ni  and  N2,  occur  within  different  compo- 
nents of  N  then  two  distinct  objects  are  merged  at  c, 
a  contradiction.  So  assume  they  are  contained  in  one 
connected  component  of  N  .  One  can  then  find  a  loop  L , 
consisting  of  cell  edgesand  cell  facediagonalsconnect- 
ing  adjacentLowson  the  boundary  faces  of  N  .that  lies 
completely  exterior  to  0  .  As  in  the  proof  of  Theorem 
2  this  loop  can  be  chosen  so  that  a  new  loop  L',  formed 
when  these  two  object  portions  merge  at  the  threshold 
c,  will  be  interlocked  with  the  loop  L .  This  shows  that 
loop  L'  isincontractibleand  thusa  changein  thefunda- 
mental  group  of  oneof  the  contours  occurs  at thesaddle 
value.  □ 
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